1 Pressure map

This chapter covers the main steps to determine the position of a bird from pressure data. This code basically implement the method presented in Nussbaumer et al.³

1.1 Read geolocator data

We read the geolocator data and crop it so that it starts on the equipment date and ends on the retrieval date.

pam <- pam_read(
  pathname = "data/0_PAM/18LX",
  crop_start = "2017-06-20",
  crop_end = "2018-05-02"
)

pam_read() currently support data from Swiss Ornithological Institute files (*.pressure, *.lux, *.acceleration) and Migrate Technology (.deg, *.lux). For other file format, please contact me (submit a Github issue or email).

1.2 Label tracks

The labeling of the track is two-fold: - Indicate the flight periods by adding a logical (TRUE or FALSE) for each acceleration measurement: pam$acceleration$ismig. - Identify outlier pressure measurement which should not be use for estimating position as a logical in pam$pressure$isoutliar.

To ensure the high level of precision needed for the pressure match, we strongly suggest to use a manually labeling through TRAINSET. The chapter labelling tracks is dedicated to this exercise and includes tips and best practices. Here, we will just walk through the main step of the workflow.

1.2.1 No acceleration data

If no acceleration data is available for your tracks, we suggest you to create a fake acceleration data, using for instance the pressure dataset:

pam$acceleration = pam$pressure
pam$acceleration$obs = 0 # you could also keep the pressure measurement for ease of identification in TRAINSET
pam$acceleration$obs[1] = 1 # trick to avoid y axis issue in TRAINSET (see https://github.com/Rafnuss/GeoPressureR/discussions/26)
pam$acceleration$ismig = FALSE

Note that while this will use the temporal resolution of the pressure measure, any resolution of activity can be used.

You won’t be able to use pam_classify(), so skip the next sub-section. You can try some other classification of migration such as change point analysis (e.g., https://kiranlda.github.io/PAMLrManual/soar.html#classify-using-a-changepoint-analysis), but in most case, it is often easier to classify directly your track manually.

1.2.2 Automatic classification of activity

We initialize the labelling file with an automatic classification of activity. We first use a k-mean clustering to group periods of low and high activity and then classify high activities lasting more than 30 minutes as migratory activities. See more possible classifications in the PALMr manual.

pam <- pam_classify(pam, min_duration = 30)

1.2.3 Edit activity on TRAINSET

Use trainset_write() to export the automatically generated classifications in a csv file, which can be opened in TRAINSET: https://trainset.geocene.com/.

trainset_write(pam, pathname = "data/1_pressure/labels/")
# browseURL("https://trainset.geocene.com/")

Print screen of the manual classification in TRAINSET. See labelling tracks for more information.

When you have finished the manual editing, export the new csv file (TRAINSET will add -labeled in the name) in /data/1_pressure/labels/) and read this file with trainset_read().

pam <- trainset_read(pam, pathname = "data/1_pressure/labels/")

1.3 Identify stationary periods

Based on the activity labelling, pam_sta() creates a table of stationary periods as illustrated below.

pam <- pam_sta(pam)
knitr::kable(head(pam$sta))

sta_id	start	end
1	2017-06-20 00:00:00	2017-08-04 19:50:00
2	2017-08-04 23:15:00	2017-08-05 19:30:00
3	2017-08-06 02:50:00	2017-08-06 19:15:00
4	2017-08-07 03:10:00	2017-08-07 19:15:00
5	2017-08-08 00:10:00	2017-08-29 18:40:00
6	2017-08-30 04:30:00	2017-08-30 18:45:00

We can visualize the pressure measurements for each grouped stationary period (symbolized by a different color). The back dots represents the pressure labeled as outlier and these data-point will not be matched.

pressure_na <- pam$pressure
pressure_na$obs[pressure_na$isoutlier | pressure_na$sta_id == 0] <- NA
p <- ggplot() +
  geom_line(data = pam$pressure, aes(x = date, y = obs), col = "grey") +
  geom_line(data = pressure_na, aes(x = date, y = obs, col = as.factor(sta_id))) +
  geom_point(data = subset(pam$pressure, isoutlier), aes(x = date, y = obs), colour = "black") +
  theme_bw() +
  scale_y_continuous(name = "Pressure (hPa)") +
  scale_colour_manual(values = rep(RColorBrewer::brewer.pal(9, "Set1"), times = 8))

ggplotly(p, dynamicTicks = T) %>%
  layout(
    showlegend = F,
    legend = list(orientation = "h", x = -0.5),
    yaxis = list(title = "Pressure [hPa]")
  )

1.4 Compute pressure maps

Now that we have clean pressure time series for each stationary period, we are ready to match each one with a weather reanalysis dataset (ERA5). To overcome the challenge of handling such a large dataset, GeoPressureR uses the API GeoPressure to perform the computation on Google Earth Engine.

Initially, it is easier and faster to query only long stationary periods (in the example below, we select only periods longer than 12hrs). You can do so by setting the pressure of the stationary periods you wish to discard to NA.

sta_id_keep <- pam$sta$sta_id[difftime(pam$sta$end, pam$sta$start, units = "hours") > 0]
pam$pressure$sta_id[!(pam$pressure$sta_id %in% sta_id_keep)] <- NA

We can now query the data on the API with geopressure_map(). A detailed description of the parameters can be found here. This will take a couple of minutes to run.

pressure_maps <- geopressure_map(
  pam$pressure,
  extent = c(50, -16, 0, 23), # coordinates of the map to request (N, W, S, E)
  scale = 2, # request on a 1/2=0.5° grid to make the code faster
  max_sample = 250, # limit the query to the first 250 data-points.
  margin = 30 # roughly equivalent to 3hPa
)

geopressure_map() returns a list of two rasters for each stationary periods. The first is the mean square error ($\textbf{MSE}$) between the pressure time series and ERA5 map. The second ($\textbf{z}_{thr}$) is the proportion of data-points in the pressure time series which correspond to an altitude that falls between the min and max altitude of each grid cell. Read more about these values and how they are computed here.

1.5 Compute probability maps

We then combine the two rasters in a single probability map using \[\textbf{P} = \exp \left(-w \frac{\textbf{MSE}}{s} \right) [\textbf{z}_{thr}>thr]\] where $s$ is the standard deviation of pressure and $thr$ is the threshold mask. Because the auto-correlation of the time series is not accounted for in this equation, we use a log-linear pooling weight $w=\log(n) - 1$, where $n$ is the number of data-points in the time series. Probability aggregation describing the influence of log-linear pooling and length of time series will be added later.

pressure_prob <- geopressure_prob_map(
  pressure_maps,
  s = 1, # standard deviation of pressure
  thr = 0.9 # threshold of the threshold proportion value acceptable
)

We use leaflet to visualize the threshold mask, mismatch map, and overall probability map for a single stationary period.

i_r <- 2
leaflet(width = "100%") %>%
  addProviderTiles(providers$Stamen.TerrainBackground) %>%
  addFullscreenControl() %>%
  addRasterImage(pressure_prob[[i_r]], opacity = 0.8, colors = "OrRd", group = "Probability") %>%
  addRasterImage(pressure_maps[[i_r]][[1]], opacity = 0.8, colors = "OrRd", group = "Mismatch") %>%
  addRasterImage(pressure_maps[[i_r]][[2]], opacity = 0.8, colors = "OrRd", group = "Threashold") %>%
  # addLegend(pal = pal, values = values(v[[i_s]][[3]]), title = "Probability") %>%
  addLayersControl(
    overlayGroups = c("Probability", "Mismatch", "Threashold"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  hideGroup(c("Mismatch", "Threashold"))

We can also visualize the probability map for all stationary periods:

li_s <- list()
l <- leaflet(width = "100%") %>%
  addProviderTiles(providers$Stamen.TerrainBackground) %>%
  addFullscreenControl()
for (i_r in seq_len(length(pressure_prob))) {
  i_s <- metadata(pressure_prob[[i_r]])$sta_id
  info <- pam$sta[pam$sta$sta_id == i_s, ]
  info_str <- paste0(i_s, " | ", format(info$start, "%d-%b %H:%M"), "->", format(info$end, "%d-%b %H:%M"))
  li_s <- append(li_s, info_str)
  l <- l %>% addRasterImage(pressure_prob[[i_r]], opacity = 0.8, colors = "OrRd", group = info_str)
}
l %>%
  addLayersControl(
    overlayGroups = li_s,
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  hideGroup(tail(li_s, length(li_s) - 1))

1.6 Compute altitude

The second operation you can perform with GeoPressureR is to compute the exact altitude of the bird $z_{gl}$ from its pressure measurement $P_{gl}$ using the barometric equation, correcting for the natural variation of pressure and temperature.

This function used is \[ z_{gl}(x)=z_{ERA5}(x) + \frac{T_{ERA5}(x)}{L_b} \left( \frac{P_{gl}}{P_{ERA5}(x)} \right) ^{\frac{RL_b}{g M}-1},\] where $z_{ERA}$, $T_{ERA}$ and $P_{ERA}$ respectively correspond to the ground level elevation, temperature at 2m and ground level pressure of ERA5, $L_b$ is the standard temperature lapse rate, $R$ is the universal gas constant, $g$ is the gravity constant and $M$ is the molar mass of air. See more information here.

To illustrate the benefit of using this equation, we will compute the bird’s altitude for its first stationary period using (1) GeoPressureR and (2) the barometric equation using standard atmosphere condition.

We first determine the position of the bird by using the most likely position using geopressure_map2path

pt <- geopressure_map2path(pressure_prob[1])

And then call the function geopressure_ts() with the subset of pressure containing sta_id==1

pressure_timeserie_1 <- geopressure_ts(pt$lon, pt$lat, pressure = subset(pam$pressure, sta_id == 1))

We can compare the altitude produced to the one computed without the correction for temperature and pressure:

Lb <- -0.0065
R <- 8.31432
g0 <- 9.80665
M <- 0.0289644
T0 <- 273.15 + 15
P0 <- 1013.25
pressure_timeserie_1$altitude_baro <- T0 / Lb * ((pressure_timeserie_1$pressure / P0)^(-R * Lb / g0 / M) - 1)

and visualize this comparison:

p <- ggplot() +
  geom_line(data = as.data.frame(pressure_timeserie_1), aes(x = date, y = altitude, col = as.factor("Corrected elevation with ERA5"))) +
  geom_line(data = as.data.frame(pressure_timeserie_1), aes(x = date, y = altitude_baro, col = as.factor("Uncorrected elevation"))) +
  labs(col = "") +
  theme_bw()

ggplotly(p) %>%
  layout(legend = list(orientation = "h", x = -0.5))

The function geopressure_ts() also returns the ground level pressure time series from ERA5 at the location specified. This is useful to check whether there is a good match between the pressure measured by the geolocator and the one at the assumed location. This operation is typically used to check the quality of the manual labelling (see labelling tracks).

1.7 Compute pressure and altitude for the path

We can repeat the computation of the pressure time series for all stationary periods. First we compute all the most likely position from the probability map of pressure.

path <- geopressure_map2path(pressure_prob)

Secondly, we can use geopressure_ts_path() which basically call geopressure_ts() in parallel for all stationary periods. We can additionally request to compute the altitude during the next flight flight with include_flight = c(0,1). Note that if a position of the path is over water, it will be moved to the closest point onshore.

pressure_timeserie <- geopressure_ts_path(path, pam$pressure, include_flight = c(0, 1))

p <- ggplot() +
  geom_line(data = do.call("rbind", pressure_timeserie), aes(x = date, y = altitude)) +
  theme_bw() +
  scale_y_continuous(name = "Altitude (m)")
ggplotly(p, dynamicTicks = T) %>% layout(showlegend = F)

col <- rep(RColorBrewer::brewer.pal(9, "Set1"), times = ceiling((nrow(pam$sta) + 1) / 9))
col <- col[1:(nrow(pam$sta) + 1)]
names(col) <- levels(factor(c(0, pam$sta$sta_id)))

p <- ggplot() +
  geom_line(data = pam$pressure, aes(x = date, y = obs), colour = "grey") +
  geom_point(data = subset(pam$pressure, isoutlier), aes(x = date, y = obs), colour = "black") +
  # geom_line(data = pressure_na, aes(x = date, y = obs, color = factor(sta_id))) +
  geom_line(data = subset(do.call("rbind", pressure_timeserie), sta_id != 0), aes(x = date, y = pressure0, col = factor(sta_id))) +
  theme_bw() +
  scale_colour_manual(values = col) +
  scale_y_continuous(name = "Pressure (hPa)")

ggplotly(p, dynamicTicks = T) %>% layout(showlegend = F)

1.8 Save

save(
  pressure_timeserie,
  pressure_prob,
  pam,
  file = "data/1_pressure/18LX_pressure_prob.Rdata"
)

Introduction

2 Light map